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This work studies the hydrodynamics of self-gravitating compressible isothermal fluids. We show 
that the hydrodynamic evolution equations are scale covariant in absence of viscosity. Then, we 
^ ] study the evolution of the time-dependent fluctuations around singular and regular isothermal 

■ spheres. We linearize the fluid equations around such stationary solutions and develop a method 
' based on the Laplace transform to analyze their dynamical stability. We find that the system is 

£\j | stable below a critical size (X ~ 9.0 in dimensionless variables) and unstable above; this criterion 

is the same as the one found for the thermodynamic stability in the canonical ensemble and it is 
associated to a center-to-border density ratio of 32.1. We prove that the value of this critical size 
is independent of the Reynolds number of the system. Furthermore, we give a detailed description 
of the series of successive dynamic instabilities that appear at larger and larger sizes following the 

■ geometric progression X n ~ 10.7", n = 1, 2, .... Then, we search for exact solutions of the hydro- 
i—l ' dynamic equations without viscosity: we provide analytic and numerical axisymmetric soliton-type 

solutions. The stability of exact solutions corresponding to a collapsing filament is studied by com- 
£NJ . puting linear fluctuations. Radial fluctuations growing faster than the background are found for all 

sizes of the system. However, a critical size (X ~ 4.5) appears, separating a weakly from a strongly 
CO unstable regime. 

i> : 
o ■ 

I. INTRODUCTION AND RESULTS 

o 

Q\ . ....... 

. Understanding the dynamics of self-gravitating fluids is a fundamental issue in astrophysics. Indeed, the structure of 
the interstellar medium as well as the formation and evolution of cosmological structures are based on it. Compressible 
sclf-gravitating fluid mechanics plays a key role in both star formation and in the formation of a fractal structure for 
the interstellar mediuml]]J. 
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Self-gravitating systems can be studied at different levels. The first is purely thermodynamic. It leads to the 
description of their states of equilibrium and the study of their stability. This approach has been followed in the case 



c/5 
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of the isothermal sphere configuration j^, ^J, ||, ||, |6| leading to the description of the gravothermal instability. This 
behaviour appears when the size of the sphere goes above a critical size, or, for a given size, when its temperature drops 
below a critical temperature. An instability develops and the system collapses to a dense core. In refs.j9|, [ll], 
k> , the statistical mechanics of the self-gravitating gas has been investigated with analytic (saddle point in the functional 
integral) as well as Monte Carlo methods. In addition to the thermodynamics, this approach provides the exact 
equation of state of the self-gravitating gas|il|. 

The dynamics can be studied from the hydrodynamic equations of a self-gravitating fluid. This is a set of non-linear 
partial differential equations and exact solutions can be found only under simplifying hypotheses. For example, the 
case of one space dimension has been studied in ref.|ll| and the axisymmetric case in ref.|Q. 

Finally, the last approach to the study of self-gravitating systems is numerical, it has been widely followed]^]. 
In this paper we study the hydrodynamics of a self-gravitating fluid. We show that the hydrodynamic evolution 
equations for a self-gravitating fluid are covariant under scale transformations in absence of viscosity. That is, if 
we scale the space and time variables by a constant factor in a solution of the fluid equations, we get again a solution 
of the equations. 

We consider the Navier-Stokes equation for an isothermal self-gravitating fluid and restrict our study to potential 
flows. Within these hypotheses we show that a Reynolds number appears as the only parameter that breaks the scale 
covariance in the Navier-Stokes equation. Moreover we find that a characteristic density enters the Reynolds number 
definition, implying that a system with no characteristic density has no characteristic Reynolds number; the flow is 
then scale covariant. 



In section III we present a general method for the analysis of linear dynamic stability of stationary solutions of 
the equations. This method is based on the Laplace transform in time of the evolution equations. The functional 
determinant of the Laplace transformed operator is studied as a function of s (the Laplace transform parameter) . Its 
zeros in the s-plane determine whether the solution is stable or not. Unstable solutions yield zeroes with positive 
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real part while purely imaginary zeros or zeroes with negative real part correspond to stable solutions (oscillating 
or damped). We apply this analysis to isothermal spheres where we first specify the initial data and the boundary 
conditions. 

In section IV , using the Laplace transform, we derive the analytic form of the fluctuations in the case of a singular 
isothermal sphere background. Then, we show that the boundary conditions build a discret spectrum of proper 
modes for the system. The stability of these modes is decided by the complex value of their pulsations in the Laplace 
transform plane. We then apply a general numerical method based on finite differences to both the singular and the 
regular isothermal sphere cases and compute the pulsations as functions of the size of the system. We show that in 
both cases the dynamic stability depends on the size of the system. Under a critical size X Cl the system is stable, 
while above it becomes unstable. In dimensionless variables this size is X C1 ~ 10.7 15 in the singular case (S is the 
short distance cutoff necessary in this case) and X C1 ~ 9.0 in the regular case. This critical size is the same as the one 
found in thermodynamic studies in the canonical ensemble (we used boundary conditions connected to the canonical 
ensemble). It is often expressed as a critical center-to-border density contrast whose value is 32.1 for the regular 



isothermal sphere. However, we also show how secondary instabilities X c 



, X Cn appear at larger sizes. These sizes 



follow the geometric progression 10.7™ (n = 1,2,.. .) in the singular case and do so asymptotically (n 3> 1) in the 
regular case. The density and velocity profiles of the two first instabilities are given. Figs. 1-3 show these profiles 
and the spectrum of modes for the regular and singular spheres. 

In section [v| we include the viscosity in the fluid equations. Then, we study the influence of the Reynolds number 
on the stability criterion. The result is that the critical sizes at which the instable modes appear are independent 
of the Reynolds number. The complicated evolution of the pulsations of the modes as functions of both the system 
size and the Reynolds number is plotted in figs. 4-6. 

In section VI we turn to time-dependent background solutions. We look for axisymmetric soliton-type solutions. 



We define a soliton variable: 



y 



/ir 



1 ± fj, c s t 



(1.1) 



where p, is an inverse characteristic length (p 1 is the Jeans' length), c s the sound speed, r the polar radius and t the 
time. We reduce the equations to a single non-linear equation and we find non-trivial solutions such as, 



±c s [ Vc 



p, r 



1 ± fJLC s t 



<f> = In 



c 



p? r 2 (l ± (i c s t) 



(1.2) 



This solution describes an expanding or collapsing filament dependin g on the choice for ±. Other background solutions 
are described. Then, explicit radial fluctuations around the solution (1.2) are computed and their stability is analyzed 
with the same method as for the isothermal spheres. Some of the modes appear to grow faster than the background 
whatever the size of the system. A specific size, X C1 ~ 4.5 5, marks the appearance of faster growing fluctuations at 
larger sizes. However, these are radial fluctuations invariant along the axis of the symmetry; consequently they have 
an infinite mass for all radial sizes of the system. Limitating the size of the system along the symmetry axis should 
provide a stab le behaviour for small enough sizes. 

Section VII is devoted to our conclusions and final remarks. 



II. THE NAVIER-STOKES EQUATION FOR SELF-GRAVITATING FLUIDS 



The evolution of a self-gravitating, isothermal fluid is described by a density field, p(r,i), a pressure field P(r,t), a 
gravitational potential U(r,t) and a velocity field v(r, t). These four fields obey four equations: the matter conserva- 
tion, the Navier-Stokes equation, the Poisson equation and the equation of state H, respectively given by 



~dt 



+ (v • V)v 



dp 
at 

VP 



+ V-(pv) 



VC/ + - V 2 v- 



W 2 U = 4ttG p, 

2 _ 9P 
Cs ~ ~dp ' 



(C 



V(V-v) , 



(2.1) 

(2.2) 
(2.3) 
(2.4) 



3 



Here C an d V stand for the kinematical viscosity coefficients and c s is the speed of sound. We assume c s to be a 
constant which corresponds to the equation of state of a perfect gas. 

^From now on, we restrict for simplicity to potential flows, which implies the existence of a velocity potential "J. 

In order to work with dimensionless variables we define n, the inverse of a characteristic length, as: 



p 2 = 



(2.5) 



where po is a characteristic density of the system. 
We then introduce dimensionless variables and fields 



(j, r 



fic s t 



$ = In - 



Po 



fi 



V x *(x, r) 



(2.6) 



Here, V x is the derivative with respect to the dimensionless coordinate x. We can take advantage of the potential 
nature of the velocity field to simplify the dissipation term in the Navier-Stokes equation. Moreover, to be able to 
eliminate the field U, we now take the divergence of Navier-Stokes equation. At the same time, to get an equivalent 
f orm ulation, we must take the rotational of Navier-Stokes equation. This last constraint yields the simple equation 
(2.9) for non-zero viscosity. It disappear for flows with zero viscosity. The system can now be described with the 
variables $ and ^ only as follows 



d T VH + (VtV)V* = -V 2 $ - e* + —V • \e^V(VH)] , 

Re 

cV* + V 2 * + V$ • W = 0, 
V$AV(V 2 $) =0, 



(2.7) 

(2.8) 
(2.9) 



JJ_ 

PoCs 



(2.10) 



Re is the analogous of the usual Reynolds number for the compressible self-gravitating flow. 

Now, let us see that the flow is scale covariant in the absence of viscosity. This can be proven as follows. Assume 
that {$(x, r), ^(x, r)} is a solution of eqs.(2.7)-(2.9). We then define the scale-transformed fields as 



$a(x,t) = $(Ax,At) +hiA 2 
*a(x,t) = A- 1 *(Ax,Ar) 



(2-11) 



Using that 



= A 



dx k 
dx k 



where yk = Xxk and f = At, we obtain that $a(x,t 



<9$ 


3$a 


dyk 


' dr 




3*A 


dyk 


' dr 


x,t), 


*a(x,t) 



- X df 
~df 



(2.12) 



transformation (2.11) generalizes to the time dependent case the scale transformation for $(x) found in rei.|T(| in the 
hydrostatic case. 



We find from eqs.(2.6) and (2.11) that the velocity field transforms as, 

v a (x,t) = v(Ax, At) . 



Un der this transformation, all terms in eq.(2.8) scale as A and in eq.(2.7) as A 2 except for the last (viscous) term in 
eq.(2.7) that scales as A. Therefore, for A — > oo scale invariance is recovered. This corresponds to the long wavelength 
limit as one could expect. 

When the flow is a fluctuation, for example, a fluctuation around an isothermal sphere with a characteristic central 
density, the scale invariance is broken. However, in this case we will see that it is recovered asymptotically at large 
distances. 
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III. THE DYNAMIC STABILITY OF ISOTHERMAL SPHERES 

We study in this sect io n th e evolution of time dependent fluctuations around an isothermal sphere. We linearize 



the fluid equations (2. 7)- (2. 9) around such exact static solutions and apply Laplace transform to solve them. 



A. Linearization 

Equation ( [2.7| ) and fl2.8p have well-known static spherically symmetric solutions: the isothemal spheres. These are 
the solutions of the equations: 

-r^ + --7 £ +e $0 =0, 3.1 

ax A x ax 

*o = , (3.2) 

where x is the dimensionless radius fir. These solutions are studied in |g, [|, |j, [l2|. Stable isothermal spheres have a 



finite radius. Solutions of eq.(B.f ) reaching infinite radial distance suppose an infinite total mass. Therefore, we shall 
consider the fluid confined in a spherical box of finite size. 

The thermodynamic stability of isothermal spheres has been extensively studied, leading to the discovery of the 
gravo- thermal instability (@, @ and @). Our aim is to investigate the dynamic instability of these states and 
compare with the thermodynamic results |Q . To achieve the stability analysis we introduce fluctuations around the 
isothermal spheres: 

$(x,r) = $ o + e0(x,r) , (3.3) 
*(x,r) = eV(x,r), (3.4) 

where e is a small parameter. 

We can then write the linearized system where only first order contributions in e are kept: 

9 T V 2 ^(x, r) = -V 2 0(x, t) - e* o(x) 0(x, r) + — V • fe"* ^ V(V 2 ^(x, r))l , (3.5) 

Re 

<9 T 0(x, t) + V 2 V(x, t) + V<J> (a;) • V^(x, r) = , (3.6) 

V<£>o A V(V 2, 0) = . (3.7) 

Since we study fluctuations around a spherically symmetric solution, it is natural to expand i/)(x, t) and </>(x, t) in 
spherical harmonics. Let us write first, 

0(x, r)=J2 <t>l,™{x, t) Yuje, <p) , (3.8) 



^(x, r) = 1>l,m(x, r) Y l>m {6, tp) . (3.9) 

l.rn 

The action of the laplacian is on </>^ m and ipi lTn is: 

V 2 ^ + ^_fcll. (3 . 10 ) 



Then, the system of equations can be written as 
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V 2 d T ^(x,T) + [V 2 + e*°]Mx,T) - -L V- [e"*"' 1 ' V(V 2 V*(^))] = , 

-Ke 

d T ct>i{x,T) + [V 2 + (d x $o) &#j(a:,T) = , 

53 v 2 ^ fl *j, m = o , 53 v 2 ^ 9^, m = o . 



(3.11) 



(3.12) 



l.m 



l.m 



We can see from the last condition that all m ^ modes must be zero. This justifies labeling the radial components 
of cj> and ip with I only. Such a cancellation ofm^O modes happens only in the linear approximation. Moreover, the 
two last conditions are trivial for I = 0; the S-wave mode is indeed decoupled from the others. This is not the case 
however for I > modes, which are coupled through the condition (3.12), if the viscosity is not zero. Consequently, 
in the viscous case, we will restrict ourselves to the study of S-wave perturbations. Since the two last conditions 
are automatically satisfied within the extend of our study, we will not mention them again. 



Analytic solutions will be given in specific cases. We analyze below the system of equations (3.11) by using Laplace 
transformation in time, showing that this is a powerful approach for a stability analysis. 



B. Laplace transform analysis 



We are interested on fields which grow slower than some exponential of time and therefore admit a Laplace transform 
for some finite s > 0. This space of solutions is rich enough for our analysis. The Laplace transforms of our fields are 
defined by: 



ipi(x,s) = 



+ OC 



dr 



+oo 



ipl(x,r) . 



M(x,s) 



The Laplace transform of the evolution equations (3.11) takes then the form 

^0,s)\ / V 2 ^(s,r = 0) \ 
$i(x,s) ) \ 4>i{x,t = Q) ) 
where 

_i_ y • ( e -*°<» VV 2 ) V 2 + e <E, °( :r ) 

V 2 + d x <I>o d x s 
Indeed, the inverse Laplace transforms can be written : 



(3.13) 



(3.14) 



M(x, s) 



s V 2 



(3.15) 



4>i{x,t) 



p 2/Kl 

e sr - 

p 2/Kl 



(3.16) 



where the contour T runs upwards parallel and to the right of the imaginary s axis. The functions 4>i(x, s) and i)i{x, s) 
in the integrand of eq.(3.16) are given by the inverse operator M.{., s) _1 acting on the r.h.s. of eq.(3.14). Hence, the 
poles of the integrand are given by the zeros Sfc, k = 1, 2, 3, . . . of the determinant of the operator A4(., s). 

Closing the contour T with a path on the Res < half-plane and computing the integral in eqs.(3.16) by residua 
we get for the fields a sum of terms associated to the poles of the integrand. Then, we can push the contour in the 
Res < half-plane towards Res — > — oo where its contribution vanishes. 



Therefore, the inverse Laplace transforms (3.16) become sum of terms depending on time through e Sk T . If S& has 
a positive real part the mode is unstable, growing exponentially with time. 



poles 



res(pi(x, s fc ) , 
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^{x,t) = e SkT resMx,s k ) . (3.17) 

poles 

The method for analyzing the stability of the system is now clear: we determine the zeroes of the determinant of our 
operator A4(., s) as a function of s for various values of I under specified boundary conditions, and we look for zeroes 
in the Re(s)> half-plane of the complex s-plane. 



C. Boundary conditions on an spherical domain 

Regularity at the origin requires, 

d x <pl( x , T )\x=o = ° . 9 x ipi{x,T)\ x=0 = . 



(3.18) 



As already mentioned, we confine our system in a spherical box to provide a large distance cut-off. We consider a 
spherical box of radius X and we require ipi an d the radial velocity to vanish on the surface of the sphere. The first 
condition has is like a gauge condition and has no physical effect, while the second prevents the fluid from exiting the 
confining box. 



-x 



= 



MX,t)=0 



(3.19) 



These conditions translate for the Laplace transforms into: 

0*01(0,*) = d x M0,s) = , d m $i(X,a)=$ l (X,a)=0 

Perturbations with non-zero total mass are in fact zero-mass perturbations around an isothermal sphere with a 
different parameter fi = They must not be included in the stability study. 

We can actually prove that the mass of the field is unchanged to first order by the perturbations. Namely, that the 
integral 



x 



(f,r) e^ x) x 2 dx 



(3.20) 



identically vanish. 

For perturbations with I ^ the integral (3. 2C) trivially vanishes upon integration over the angles. For I = 0, we 
can rewrite the Laplace transform of eq.(3.6) as follows 



s fa(x, s) + e -*° (x) V [e* o(3;) V^o(a;, a) 



Hence, integrating over the sphere, 



X M^s) e*°^ x 2 dx = ~e^ ^° (X ' S) 



dx 







x=X 



= 



x=0 



where we used eqs.( |3.1£| )-( |3.1S ). 



I n ord er to study the solutions of eqs.(3.14) with the boundary conditions eq.(3.1S)-(3.19), let us first reformulate 
eq.( 3.14 ). We will restrict this study to I = 0. 

Introducing f(x,s) 



^f 2 - we can reduce eq. (3.14) to: 



d_ 

dx 



-*oO) 



Re 



■ d 2 f 




dx 2 


X v 



1+ R~e e 



-*o(») 



1 - 

2 d<P 
x dx 



-*d(x) 



Re 



dx 
dx 2 



dx 



0) 



^-(x,s) + ^-f(x,s) 
dx dx 



(3.21) 
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We can integrate the total derivative and drop the integration constant thanks to eqs.(3.1£)-(3.19). We find, 



1 + —e-foix) 
Re 



dx 2 



±(i + ^ e -*o(*n + ^o 



,-*o(a0 



Re 



Re 
2 

x dx 



dx 



dx 



f(x,s) 



A X , T = )=0 
dx 



The boundary conditions for f(x,s) follow from eqs.( [3.18 )-( [3.1S| ): 

f(0,s) = f(X, s) = 



Let us first analyze this linear problem in the zero viscosity case, namely Re — oo. Eq.(3.22) thus becomes, 





dx 2 



"2 d<I> ~ 




x dx 


dx 



X 



2 2 d$ 

2 x dx 



f(x,s) + S ~l^-( X : T = 0) 



(3.22) 



(3.23) 



(3.24) 



This ordinary second order differential equation takes a Schrodinger-type form upon the transformation 



e -i * (a:) 
f(x, s) = - w(x, s) 



We find for the homogeneous part 



d 2 w 
dx 



u, w . . 

—( x > s > + 



1 ( d$ 

dx 



2 d$ Q 1 

x dx 2 



•I'D 



w(x,s) — —s 2 w(x,s) 



(3.25) 



Here. 



V( \ - 2 JL 1 



2 d£o _ l e$0 
x dx 2 
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/d£o _ 


2 X 


\~\ 






2 


V c?x 


x . 




\ dx . 





plays the role of the 'potential' and — s 2 the role of the 'energy' eigenvalue in the Schrodinger-type equation. 
Notice that the potential is singular and repulsive at short distances since 

V(x) x =° ^+0(1) , 

while it is attractive for distances where ^ tt ~ —2/x. In the case of a potential attractive at short distances, or 
in the case of an infinite domain, a continuous part could appear in the spectrum. In our case (finite domain and 
repulsive potential), the spectrum will be discrete (see ref. |l7j p 445 for a mathematical study). Notice however that 
the numerical methods used in the sequel can handle both continuous and discrete sprectrum. 

Thus, imposing the boundary condition ui(0, s) = w(X,s) = yields a discrete spectrum of eigenfunctions w n (x) 
wit h neg ative eigenvalues s 2 t plus possibly some eigenvalues with positive corresponding to the 'bound states' in 
eq.( 3.25| ). This ensemble of eigenfunctions form a comp lete set. 

We can now solve the inhomogeneous equation (3.24) by expanding in the eigenfunctions 



fn(x) 



*o(a0 



w n (x) 



of the homogeneous problem (3.25) with boundary conditions (3.23) and normalized as 



x 



We thus write 



x 2 e*°W f n (x) f k (x) dx = S nk 
f(x,s) =^c n (s) f n (x) 



(3.26) 



where the c n (s) are for the moment arbitrary. Inserting eq.( p.2q ) into eq.(3.24), multiplying by fk{x) and integrating 
over x yields, 



f(x,s) = s V / 

: Jo 



x x' 2 e*^') dx! dj) 

— 2 fn(x) fn{X ) — (x , T = 0) 

s z — st dx' 
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Thus, f(x,s) exhibits poles in s 2 at each eigenvalue s 2 . Positive s 2 describe instabilities as discussed at the end of 
sec. IIIC since such poles imply exponentially growing fluctuations as e' s,l ' T . Negative eigenvalues s 2 n = — fc 2 yield 
oscillating behaviour e lknT . 



In the viscous case eq.( 3.22 ) can be reduced to a 'Schrodinger type' equation with an energy dependent potential 
upon the transformation 



f(x,s) 



e 2 



5*0 (x) 



w(x, s) 



where 



We then find for the homogeneous part, 



(3.27) 



dx 



2 {X, s) + 



2 d$ 



\ dx ) 

Aa 2 (x) x dx 2 a s (x) Re 2 a 2 (x) 



' d<t> 
, dx 



-#q(x) 



w(x, s) 



a 2 s {x) 



w(x, s) 



A more mathematical study of the fluctuation equations can be done here in the spirit of ref. . 



IV. THE NON- VISCOUS ISOTHERMAL SPHERES 



A. Analytical study 



The singular isothermal sphere is an analytic solution to eq (3.1) given by 



$o = In -j 

x z 



(4.1) 



This solution has an infinite total mass and produces infinite pressure at the center. To avoid this divergence we 
introduce a short distance cutoff 5. The spherical box provides the long distance cutoff. We study the singular sphere 
background as an analytic model for the regular sphere. In this spirit we rewrite the boundary conditions as: 



d x ipi(x,s) 
d x ipi(x,s) 



x—8 



x=X 







d x $l(x,s) 



X — 5 

, $ l (X,s)=0. 



(4.2) 



Conditions at x — S are straightforward translations of conditions at x = 0. 

The advantage of this a nalyt ic background solution is that we can write down analytic expressions for the fluctua- 
tions. Indeed for / = eq.( 3.25 ) reduces to: 



dx 2 



f(x,s) = 0, 



(4.3) 



where f(x, s) = w(x, s) in this case, is the dimensionless radial velocity of the / = mode. In dimensionful variables, 
the radial velocity fluctuations for I = are: 



v r (r,t) = c s (A ie lut + A 2 e~ wt ) Vkr [fli J v (kr) + B 2 J- V (kr)] , 



= ik 



= k 2 



iV7 



(4.4) 



where J v (z) are Bessel functions. The corresponding fluctuation cj>(r,t) takes the form 
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,t) = (A' l( 



A'e~ lult 



) Vk7 



D 



dr 



■(kr) 



(4.5) 



The boundary conditions (3.23) yield a system of two homogeneous equations for the constants B%, B 2 , 

B 1 J v {kX) + B 2 J- u (kX) = 

B x J u {k5) + B 2 J_ v (k6)=0 

Therefore, in order to have a non-trivial solution we impose the trascendental equation, 

J v {kX) _ J u {k5) 
J- v {kX) ~ J_„(fc<5) 



(4.6) 



(4.7) 



This equation is fulfilled f or sp ecific values of fc, yielding a discret spectrum for the modes. 

Real solutions k n of eq.(^7) can be obtained asymptotically using the Mac Mahon expansion |l5|. Depending on 
the value of X/5 there can be solutions with real or with imaginary k n . The threshold of instabilities precisely appears 
when the eigenvalues go through k = turning from real to imaginary. In such limit eq.(4.7) yields 



V7. X c 
— log — — 

2 B 5 



with the exact solutions 



X r , 2jlL 

-p- = e^ = (10.749...)' 







I = 1,2,3,. 



Each time 4- crosses each of these thresholds upwards, a new instable eigenvalue appears. The first instability occurs 



for 1 = 1. These points exactly coincide with those found in ref . 1 1 1 



We have c heck ed that the general method that we will describe in section (IVB) gives numerically the same results 
as eqs.(0)-(0). 



B. Numerical method 



The general numerical method, consists in replacing the derivatives by finite differences turning eq.(3.14) and the 
boundary conditions into a N x N (N ^> 1) linear system. More specifically, <j> and if) as functions of x are represented 
by their values on a discret set of points {xk, 1 < k < N}. The differential operator are then represented by finite 
differences. Evaluating the differential equation in each point produces an algebraic equation. Boundary conditions 
provide the missing equ ation s at the boundary of the domain, where the differential operator cannot be evaluated. 
As explained in section III B , we then evaluate the determinant of the system in a given domain of X/5, s and Re 
and look for cancellations. 

As the parameter domain is huge (four-dimensional) , we must limit the number N in the determinant. It has been 
set to 40 in most computations. However, different values (up to 200) have been used to investigated the stability 
of the results around randomly chosen cancellations of the determinant, and more specifically in the region of the 
stability-instability transition. No ill-behaviour has been detected. 



C. Stability of the singular sphere 

Cancellations of the determinant of the differential operator happen for discret values of the parameter s. If one of 
these eigenvalues s has a positive real part, then the associated mode grows with time and the system is unstable. It 
is known from thermodynamic studies that for a given central density, the stability depends on the size of the sphere. 
Thus, we will investigate for the zeroes of the determinant in a three-dimensional domain (s, X) eCxi?. The zeroes 
are shown on fig.|l| for 1 = 0. 

The first result of fig|l| is that it shows a symmetry around the Re(s)= axis; we have a set of coupled modes 
with opposite pulsations. Then, we have checked that the pulsations of the spectrum of modes have purely imaginary 
values for X/5 smaller than a critical value X Cl /5. This range of X/5 values in mostly outside the domain pictured on 
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since the behaviour of the system is simple there (however \n(X Cl /S + 1) ~ 2.4 is visible on the figures). In this 
regime the modes are oscillating stationary modes and the system is stable. As X increases above X Cl , the size of the 
system grows and a first mode encounters a branching point and develops a real positive pulsation. This means that 
the mode grows exponentially and the system becomes unstable. Then, for successive critical values X Cil i = 1,2, . . . 
corresponding to larger and larger sizes of the system, new modes become unstable. It is interesting to mention that 

the critical values appear when the rate of the larger and lower cutoffs is X Cn /5 ~ 10.7", n = 1,2, This geometric 

progression has been found in the thermodynamic study in ref. [0. When X goes to infinity, an infinite number of 
unstable modes appear, which is not surprising since the system then has an infinite mass. Obviously, this limit is 
unphysical. 

Non-radial perturbations have been investigated (1^0). These perturbations are always stable and oscillating. 



D. Stability of the regular isothermal spheres 



Regular solutions of eq.(3T) have been widely studied in spite of the fact that they do not have an analytic 
expression. In this case it is not necessary to introduce a short distance cutoff in the study of the fluctuations since 
the background solution is regular at the origin. We keep considering the solution in a large but finite sphere of radius 
X. Thus we will use the same boundary conditions as in the singular case, taking (5 = 0. 

We again apply the general numerical method described above. We compute the spectrum of pulsations of the 
modes for varying values of the size X of the sphere (see fig.||). 

The first point is that the behaviour is very similar to the singular case. The interpretation of figj2| is the same as 
that of fig.0. This means that the singular isothermal sphere with short distance cutoff is a good analytic analogue of 
the regular isothermal sphere. The only difference in the behaviour of the fluctuations is in the specific values of the 
X Ci . In the regular case, the instability appears at a smaller size (X Cl /5 ~ 9.) than in the singular case (X Cl /S ~ 10.7). 
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FIG. 2: Spectrum of the modes around the regular isothermal spheres. Explanations are the same as in Fig hi 



As expected, the value X C1 /S ~ 9.0 corresponds to a ratio of 32.1 between the central and border densities of the 
background. This critical ratio is the same as in the study of the thermodynamic stability in the canonical ensemble^. 
Finally the numerical study suggests that the ratio 10.7 between two successive X Ci is asymptotically obtained at 
large X in the regular case, as it must be. 

E. Profile of the instable modes 

Having established that the dynamic instability appears for the same center-to-border density contrast as in the 
thermodynamic instability, we now check if the profile of the first instable mode matches the profile of the thermody- 
namic instability, and we investigate the profile of other instable modes. 

Those profiles arc given, when we use the finite differences method, by the coordinates of the vector of the one- 
dimensional kernel of our operator (boundary conditions included). For a given size, we are thus able to determine 
the profile of each mode associated with a zero of the determinant for specific values of s. On fig. |^ the profiles are 
given for the first and second instabilities for sizes just above their critical sizes. As can be seen, the first instability is 
a spherical collapse forming a dense core. Its density profile is similar to the profile of the thermodynamic instability 
. The second instability combines a core collapse and the ejection of a shell. 

V. EFFECT OF THE VISCOSITY ON THE DYNAMICS AND STABILITY CRITERIA 



It is now interesting to establish how viscosity alters the dynamics of the fluctuations. We first proceed to an 
analytic study of the fluctuations and then perform a thorough numerical analysis. 
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FIG. 3: Radial profiles of the first two instabilities, at sizes just above their emergence, plotted for arbitrary amplitudes. 



A. Analytic study 

Here we will study the fluctuations around the singular isothermal sphere with a non-zero viscosity. We consider 
only radial perturbations since it simplifies the equations and contains the essential information about stability. Using 
the variables: 

y = sx, f(y) = Wip{sx), (5.1) 



where / is the (dimensionless) radial velocity, the system (3.14) becomes 



s Re 2 J dy 2 s Re dy \y 2 sRe J 



f = . (5.2) 



It is interesting to notice that all the dependence on the parameters is through the combination s Re, which is the 
turbulent-like Reynolds number associated to a mode of pulsation s. The asymptotic behaviour of the solutions at 
large distance is: 



y > 1 , y > s Re : f(y) = y * 



(5.3) 
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For sRe > — 1, the component with the upper sign grows at large radius while the other component decays. For 
negative sRe < — 1, both solutions decay at large radius. The interval — 1 > s Re > — 9/8 is the only range giving 
strictly decaying behaviours at large radius. In the other cases (s Re < —9/8 or s complex) we have spatially oscillating 
and decaying behaviours. It is easy to check that at small distances one recovers the zero- viscosity solutions. This, of 
course, does not take the boundary conditions into account. We provide in the next section quantitative results from 
numerical calculations. 



B. Numerical study 



The introduction of the viscosity, although formally cumbersome, fits with no particular difficulty in the finite 
difference scheme . In particular, it does not change the actual order of the system of differential equations as can 
be seen from eq.( 3.21 ). The spectrum of modes for a continuous range of the larger cutoff X can be computed for 
different values of the Reynolds number. These results are shown on fig.|4| and fig.[| We show in fig.|6| the evolution of 
the spectrum at fixed sizes, for a continuous range of values of the Reynolds number. 

Obviously, the diagrams are more complex in the viscous case than in the ideal case, especially as we go to Reynolds 
numbers close to one. Reynolds numbers as small as 0.1 have been investigated but are not plotted since they show no 
specificity and are unlikely values for a self-gravitating gas. We have three main kind of modes: oscillating-decaying 
modes, decaying modes, and growing modes. However, all the values of the critical sizes associated with the onset of 
instabilities are independent of the Reynolds number. This can be guessed on fig.||, and checked on fig ^ for X Cl . 
The stability of the isothermal sphere does not depend on the viscosity. Besides, figs^ and ||| clearly show that the 
mode associated with the first instability is much less affected by the viscosity than the others. 
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The real parts of the modes are plotted on fig.|6| as functions of the Reynolds number with logarithmic scales, for 
two different sizes: 10.5 and 11.6 . In the first diagram all modes are stable independently of the Reynolds number 
and in the second there is exactly one unstable mode for all values of the Reynolds number. This is a confirmation 
that the dynamic stability criterion depends only on the size of the system and not on the value of the Reynolds 
number. 



VI. SOLITON-TYPE SOLUTIONS IN THE SELF-GRAVITATING GAS 

The appearance of solitons in an isothermal self-gravitating gas has been studied in one spatial dimension in ref. |l3j . 
In this case, the system could be reduced to a sine-Gordon equation. Here we will study the case of axisymmetric 
soliton-type solutions in 3 + 1 dimensions. 

A. Soliton-type equations 



We will consider z-independent solutions of the evolution equations (2.7 )-( 2.1C ) without viscosity (Re = oo). It is 
useful to make the following change of variables, 



M) 



ro(t) 



(6.1) 



where tq (i) is a time dependent characteristic length. The simplest non-trivial choice for ro (t) is a linear dependence 
on time. Using the characteristic parameters of the system this leads the soliton variable: 



V 

We write the potential and radial velocity field as: 

m 



fj, r 



1 ± fj, c s t l±r 



$ = In 



Then, the matter conservation equation ( |2.§| ) yields 

h'(y) 1 1 1 m =c 

Kv) y f(y) 

which can be immeadiately integrated as 

f(y) h(y) y = TC 



city 



(6.2) 



(6.3) 



(6.4) 



The symbol =p in front of C is introduced for later convenience. Each possibility is associa ted with o ne o f the choices 
for the soliton variable. C going to zero means that the density goes to zero. Inserting eqs^Eyi) and (6.4) in the Eulcr 
equation of motion (EJ) gives the non linear equation: 



h(y(hh'±h- 



= ±C 



(6.5) 



It is possible to derive some analytic solutions of this equation. 



B. Analytic solutions 



The simplest solution to eq.f^) is h(y) = fvC In terms of physical fields this produces non-trivial solutions: 



*o(r, t) =±c s 



C 



/i 1 ± c s t 



, $ (r,i) = In 



C 



/j?r 2 (l ± /j.c s t) 



In- 



In(VCy) 



(6.6) 
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These two solutions are filaments, one expanding and the other collapsing in a finite time. This solution is invariant 
under the following scale transformations 

C^AC, c s ^ , /x — ► A /i . (6.7) 
A 



Hence, in the calculation of the fluctuations Jsec. VI C], we can set for example \JC = 1. 

Guided by numerical integration (s ee fi g. M), another type of solution can be found for low densities. It is possible 
to build a perturbative solution of eq.( |6.5| ) around the exact solution C = 0, h{y) = ^fy. Since the system is ill-defined 
at C = 0, we will consider C<1 and we will develop h(y) perturbatively in C. 

h(y) = T y[l + Ch 1 (y) + 0(C 2 )] . (6.8) 

Then h\(y) obeys the simple equation: 

y 2 h 1 + y(y 2 -l)h' 1 =±ln- , (6.9) 

y 

which has the general solution, 




hl{y) = 7W=f l A± 1 : * ] ' (,Kl01 

We see that in general h\(y) is singular at y = ±1, that is, on the wavefront x = |1 ± t\. Regular solutions can be 
obtained by appropriate choices of A and B. We plot in fig.(0) a regular solution hi(y) for B = 1. 
We obtain from eqs.(|6.8D and (6.10) the physical fields: 

V = ±c s Cyh 1 (y) , $ = ln(-^[l-(7/i 1 (y)]l . (6.11) 



It is important to notice that this perturbative development fails at short distances (y — > 0) where the density is not 
small anymore. 



Fluctuations around a soliton-like solution 



We intend to assess the stability of the solution (6J3). By stability we mean that the growth of the fluctuations 
should be slower than the growth of the background. 

We compute here the explicit form of the radial fluctuations around the soliton-like solution described by eq.( |6.6| ). 
The collapsing filament case is studied. We define the two independent variables: 



y = , z = 1 - t . (6.12) 

1 — T 



We consider the perturbative expansion 



city dtyn 

${x,T) = $ + ef(y,z), ( x ,T) = —^-c s eg{y,z). (6.13) 

ax ax 



We can then linearize the evolution equations and we find 



VCy^-y Z fl--(-Vc + y) d /- 2 d /- 9 = y^l + ^ + VCf, (6.14) 
oy z oyoz ay oz ay 1 ay 



n °f , 9f dg 

C— + z— + — = 0. (6.15) 
ay oz ay 

These equations are homogeneous in z so their solutions can be written as a sum of solutions of the type: 



f{y,z) = f(y) z a , g(y,z) = g{y) z a , a £ C 



(6.16) 



16 



Inserting these expressions in eqs.(6.14) and (6.15) and eliminating g, gives an equation for /: 



x(C-l)f"+ 2(C-l)-(2a + l)VCx f" + -(4a + 3)VC + a(a + l)x f + 2(a + l)a f = 



(6.17) 



This equation can be integrated in the general case. However, according to the transformation (6.7) we can choose 
the case \fC = 1. In this case the solutions take the following form: 



f(y) = B 1 $ (Ax, vx,Py) + B 2 yi® (A 2 , u 2 ,/3y) 



where, 



Ai = 2 , v x =2 



1 



l + 2a 



A 2 = 



2a 



2a + 1 



V 2 



2a + 1 



(6.18) 



(6.19) 



a(a + 1) 
2a + 1 



2a + 1 



(6.20) 



Here $(A, f, y) stands for the confluent hypergeometric function. It is regular at y — and grows exponentially for 
y — > oo. Expression (6.18) is degenerate for a = — |: 



a 



(6.21) 



According to ( 6.16| ), it diverges when the time r reaches 1. 
Another degenerate case is a = — 1, 



a = -1 



f(y) = V 



1 



+ S 2 



(6.22) 



Interestingly, i?e(a) 6 [— 1, — i] are the only values of a where the fluctuations are regular at y = 0. But again the 
fluctuations diverge for r — ► 1. 

Since the fluctuations are ill-behaved for y — ► oo and often at y = 0, we resort to the same method as in the singular 
isothermal sphere. We define a large distance cutoff X and a small distance cutoff 5 for y and we set the fluctuation 
to zero on these walls. It should be noted that these walls move with the soliton, and that they cannot be considered 
as physical since the constant component of the background velocity field flows through the walls. However, those 
cutoffs allow us to search for fluctuations as eigenmodes in the same way we did for isothermal spheres. These modes 
appear for special values of a defined here again by the cancellation of a determinant: 



$(\i,lii,/38)X'i§(\ 2 ,ii2,0X)-9(\uiii,l3X)F$(\2,na,l3S) = (6.23) 

These cancellations must be investigated in the complex plane for a and for various values of The results appear 
on fig. ||. The main result is that modes with Re(a) < exist for all values of ^ investigated. These modes have a 
faster growth than the background. Besides, 4- ~ 4.5 is a special value for the radial size of the system. Below this 
value all the unstable modes have the following time behaviour 

(l-r) Ql cos[a 2 In(l-r) + 0] 

that is, growing oscillations. Above this value 4- ~ 4.5 two modes appear that have a pure power law behaviour in 
time. These modes have larger exponents than the others, that is to say they grow faster. Moreover, just above this 
critical value, the exponent in the power law seems to go to —oo, which indicates a strong instability. 

We must keep in mind the fact that this is a partial analysis. Indeed our fluctuations are invariant along the z 
symmetry axis. This means that they have an infinite mass. If the size of the system along the z symmetry axis 
was bounded and fluctuations with r and z dependence computed, a stable region would certainly appear for small 
enough sizes. We believe that the value 4 ~ 4.5 is actually connected to the onset of a radial instability, while below 
this value the growth of the modes is due to the axial infinite size of the fluctuations. 
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VII. CONCLUSIONS 



In the first part of this work we have described a general method to study the dynamic stability of stationary 
configurations of self-gravitating fluids and applied it to the isothermal spheres. The interest of this method is that it 
can be applied to any stationary field. Moreover, we have shown that it is perfectly able to handle viscosity. We have 
applied this dynamical stability analysis to a stationary solution, the isothermal sphere. In this case we have found 
a series of sizes (X Cn ~ 10.7™, n = 1, 2, . . .) associated with the appearance of instable modes. The size associated to 
the first instable mode matches the critical size found in the thermodynamic theory. Moreover, we have shown that 
the values of these sizes are independent of the Reynolds number. 

In the second part of this paper we studied exact dynamical solutions. We have presented elements of a soliton 
theory for a self-gravitating perfect isothermal fluid with axial symmetry. This method provides new dynamic exact 
solutions. We have analyzed the stability of one of these solutions, a collapsing filament, with a method similar to 
the method used for the isothermal sphere. Here again we have found that a critical radial size appears to define two 
regimes. Above this size the system is unstable. Below, it is weakly unstable and may in fact be stable if the axial 
size of the system was bounded. 

These two studies are complementary descriptions of the dynamics of the self-gravitating fluid. One description 
deals with a stationary background and the other with a dynamical one. Although the first description obeys spherical 
symmetry and the second obeys axial symmetry, they are closely connected by their structure. Indeed, the density 
field of eq.(|6.6|) is the sum of two terms: ln(2//i 2 r 2 ) (identical to the stationary spherical solution) and a second term, 



In 



function of the soliton-time variable alone. In the stationary case, the instability appears when a mode 



is growing; in the dynamical case, it appears when a mode is growing faster (or decaying slower) than the evolving 
background solution. In this respect, the stability of the solution in both cases is governed by Jeans- like instabilities 
whose emergence depends only on the size of the system. We showed that the existence of a critical size does not only 
apply to stationary solutions but to dynamic solutions as well. 

Finally, we would like to add that it is probably possible to give a more general formulation of the soliton-type 
equations of the system, following on the results in ref. |l3| in the one dimensional-case. For example, the choice of 
the soliton variable in our work, eq.([T^), is probably a particular case of a more general change in the variables. In 
any case, the soliton methods provide a powerful approach to self-gravitating fluid dynamics which certainly deserves 
more investigation. 
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FIG. 4: Spectrum of viscous modes around the singular isothermal sphere. The real and imaginary parts of the pulsations of 
the modes are plotted in 3d-views for a range of values of the size X/5, and four values of the Reynolds number. All scales are 
logarithmic. 
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FIG. 5: Upper views of fig. |4j The real parts of the pulsations of the modes appear as functions of the ratio X/S. Positive real 
parts appear for the same values of X/8 (ln(l + X/S) ~ 2.4), independently of the Reynolds number. 
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FIG. 6: The real parts of the pulsations of the modes are plotted as functions of the logarithm of the Reynolds numbers for 
fixed sizes: 10.5 on the left and 11.6 on the right. 




21 




